﻿using System;
using System.Collections.Generic;
using System.Text;

namespace MLForgeSharp.Models.UnsupervisedLearningModels.DimReduction
{
    /// <summary>
    /// 主成分分析
    /// </summary>
    public class PCAModel
    {
        private double[,] components; // 主成分矩阵
        private double[] mean; // 数据的均值

        public PCAModel()
        {
            components = null;
            mean = null;
        }

        // 计算数据的均值
        private double[] ComputeMean(double[][] data)
        {
            int numSamples = data.Length;
            int numFeatures = data[0].Length;
            double[] mean = new double[numFeatures];

            for (int i = 0; i < numSamples; i++)
            {
                for (int j = 0; j < numFeatures; j++)
                {
                    mean[j] += data[i][j];
                }
            }

            for (int j = 0; j < numFeatures; j++)
            {
                mean[j] /= numSamples;
            }

            return mean;
        }

        // 计算协方差矩阵
        private double[,] ComputeCovarianceMatrix(double[][] data, double[] mean)
        {
            int numSamples = data.Length;
            int numFeatures = data[0].Length;
            double[,] covarianceMatrix = new double[numFeatures, numFeatures];

            for (int i = 0; i < numFeatures; i++)
            {
                for (int j = 0; j < numFeatures; j++)
                {
                    double sum = 0.0;
                    for (int k = 0; k < numSamples; k++)
                    {
                        sum += (data[k][i] - mean[i]) * (data[k][j] - mean[j]);
                    }
                    covarianceMatrix[i, j] = sum / (numSamples - 1);
                }
            }

            return covarianceMatrix;
        }

        // 计算特征向量
        private double[,] ComputeEigenvectors(double[,] matrix)
        {
            int numFeatures = matrix.GetLength(0);
            double[,] eigenvectors = new double[numFeatures, numFeatures];

            // 这里使用简单的幂迭代法计算特征向量
            Random rand = new Random();
            for (int i = 0; i < numFeatures; i++)
            {
                double[] vector = new double[numFeatures];
                for (int j = 0; j < numFeatures; j++)
                {
                    vector[j] = rand.NextDouble();
                }

                // 幂迭代
                for (int iter = 0; iter < 100; iter++)
                {
                    double[] newVector = new double[numFeatures];
                    for (int j = 0; j < numFeatures; j++)
                    {
                        for (int k = 0; k < numFeatures; k++)
                        {
                            newVector[j] += matrix[j, k] * vector[k];
                        }
                    }

                    // 归一化
                    double norm = 0.0;
                    for (int j = 0; j < numFeatures; j++)
                    {
                        norm += newVector[j] * newVector[j];
                    }
                    norm = Math.Sqrt(norm);

                    for (int j = 0; j < numFeatures; j++)
                    {
                        vector[j] = newVector[j] / norm;
                    }
                }

                // 保存特征向量
                for (int j = 0; j < numFeatures; j++)
                {
                    eigenvectors[j, i] = vector[j];
                }
            }

            return eigenvectors;
        }

        // 训练模型
        public void Train(double[][] data, int numComponents)
        {
            // 计算均值
            mean = ComputeMean(data);

            // 中心化数据
            double[][] centeredData = new double[data.Length][];
            for (int i = 0; i < data.Length; i++)
            {
                centeredData[i] = new double[data[i].Length];
                for (int j = 0; j < data[i].Length; j++)
                {
                    centeredData[i][j] = data[i][j] - mean[j];
                }
            }

            // 计算协方差矩阵
            double[,] covarianceMatrix = ComputeCovarianceMatrix(centeredData, mean);

            // 计算特征向量
            double[,] eigenvectors = ComputeEigenvectors(covarianceMatrix);

            // 选择前 numComponents 个主成分
            components = new double[numComponents, eigenvectors.GetLength(1)];
            for (int i = 0; i < numComponents; i++)
            {
                for (int j = 0; j < eigenvectors.GetLength(1); j++)
                {
                    components[i, j] = eigenvectors[i, j];
                }
            }
        }

        // 数据降维
        public double[][] Transform(double[][] data)
        {
            int numSamples = data.Length;
            int numComponents = components.GetLength(0);
            double[][] transformedData = new double[numSamples][];

            for (int i = 0; i < numSamples; i++)
            {
                transformedData[i] = new double[numComponents];
                for (int j = 0; j < numComponents; j++)
                {
                    double sum = 0.0;
                    for (int k = 0; k < data[i].Length; k++)
                    {
                        sum += (data[i][k] - mean[k]) * components[j, k];
                    }
                    transformedData[i][j] = sum;
                }
            }

            return transformedData;
        }
    }

    // 示例程序
    public class PCAModelExample
    {
        public PCAModelExample()
        {
            // 示例数据
            double[][] data = new double[][]
        {
            new double[] { 1.0, 2.0, 3.0 },
            new double[] { 4.0, 5.0, 6.0 },
            new double[] { 7.0, 8.0, 9.0 },
            new double[] { 10.0, 11.0, 12.0 }
        };

            // 创建模型
            PCAModel model = new PCAModel();

            // 训练模型
            model.Train(data, numComponents: 2);

            // 数据降维
            double[][] transformedData = model.Transform(data);

            // 打印降维后的数据
            System.Console.WriteLine("降维后的数据:");
            for (int i = 0; i < transformedData.Length; i++)
            {
                for (int j = 0; j < transformedData[i].Length; j++)
                {
                    System.Console.Write(transformedData[i][j].ToString("F4") + " ");
                }
                System.Console.WriteLine();
            }
        }
    }
}
